A numerical analysis of finite Debye-length eff'ects in induced-charge electro-osmosis 
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For a microchamber filled with a binary electrolyte and containing a fiat un-biased center electrode 
at one wall, we employ three numerical models to study the strength of the resulting induced- 
charge electro-osmotic (ICEO) fiow rolls: (i) a full nonlinear continuum model resolving the double 
layer, {ii) a linear slip-velocity model not resolving the double layer and without tangential charge 
transport inside this layer, and {Hi) a nonlinear slip-velocity model extending the linear model 
by including the tangential charge transport inside the double layer. We show that compared to 
the full model, the slip-velocity models significantly overestimate the ICEO fiow. This provides a 
partial explanation of the quantitative discrepancy between observed and calculated ICEO velocities 
reported in the literature. The discrepancy increases significantly for increasing Debye length relative 
to the electrode size, i.e. for nanofiuidic systems. However, even for electrode dimensions in the 
micrometer range, the discrepancies in velocity due to the finite Debye length can be more than 
10% for an electrode of zero height and more than 100% for electrode heights comparable to the 
Debye length. 

PACS numbers: 47.57.jd, 47.61.-k, 47.11. Fg 



I. INTRODUCTION 

Within the last decade the interest in electroki- 
netic phenomena in general and induced-charge electro- 
osmosis (ICEO) in particular has increased significantly 
as the field of lab-on-a-chip technology has developed. 
Previously, the research in ICEO has primarily been con- 
ducted in the context of colloids, where experimental and 
theoretical studies have been carried out on the elec- 
tric double layer and induced dipole moments around 
spheres in electric fields, as reviewed by Dukhin ^ and 
Murtsovkin 0. In microfluidic systems, electrokineti- 
cally driven fluid motion has been used for fluid manip- 
ulation, e.g. mixing and pumping. From a microfab- 
rication perspective planar electrodes are easy to fabri- 
cate and relatively easy to integrate in existing systems. 
For this reason much research has been focused on the 
motion of fluids above planar electrodes. AC electroki- 
netic micropumps based on AC electroosmosis (ACEO) 
have been thoroughly investigated as a possible pumping 
and mixing device. Experimental observations and the- 
oretical models were initially reported around year 2000 
[3, 4, 5, 6], and further investigations and theoretical ex- 
tensions of the models have been published by numerous 
groups since 0, i, S 0, [U El ■ Recently, ICEO flows 
around inert, polarizable objects have been observed and 
investigated theoretically [ll [11 [H, [H, [13, [ill ■ For a 
thorough historical review of research leadin g up to these 
results, we refer the reader to Squires et al. [IjI and ref- 
erences therein. 

In spite of the growing interest in the literature not 
all aspects of the flow-generating mechanisms have been 
explained so far. While qualitative agreement is seen be- 
tween theory and experiment, quantitative agreement is 



often lacking as reported by Gregersen et al. [ll| , Harnett 
et al. [ll], and Soni et al. [l^. In the present work we 
seek to illuminate some of the possible reasons underlying 
these observed discrepancies. 

ICEO flow is generated when an external electric field 
polarizes an object in an electrolytic solution. Counter 
ions in the electrolyte screen out the induced dipole, 
having a potential difference C, relative to the bulk elec- 
trolyte, by forming an electric double layer of width Ad 
at the surface of the object. The ions in the diffuse part 
of the double layer then electromigrate in the external 
electric field and drag the entire liquid by viscous forces. 
At the outer surface of the double layer a resulting effec- 
tive slip velocity VgUp is thus established. Many numer- 
ical models of ICEO problems exploit this characteristic 
by applying the so-called Helmholtz-Smoluchowski slip 
condition on the velocity field at the electrode surface 
[20! . [2l| . Generally, the slip-condition based model re- 
mains valid as long as 
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where k'B,T/{Ze) is the thermal voltage and Oc denotes 
the radius of curvature of the surface [H. The slip- 
velocity condition may be applied when the double layer 
is infinitely thin compared to the geometrical length scale 
of the object, however, for planar electrodes, condition 
is not well defined. In the present work we investi- 
gate to what extent the slip condition remains valid. 

Squires et al. 13] have presented an analytical solu- 
tion to the ICEO flow problem around a metallic cylin- 
der with radius Oc using a linear slip-velocity model in 
the two dimensional plane perpendicular to the cylin- 
der axis. In this model with its infinitely thin dou- 
ble layer, the surrounding electrolyte is charge neutral, 
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and hence the strength of the ICEO flow can be defined 
solely in terms of the hydrodynamic stress tensor cr, as 
the mechanical power Pmoch = !f\r\=a ''^ ' ' 'I'siipda ex- 
erted on the electrolyte by the tangential slip-velocity 
■Wsiip = Ucot, where h and t is the normal and tangential 
vector to the cylinder surface, respectively. In steady 
flow, this power is equal to the total kinetic energy dissi- 
pation Fkin — <|r|(^«^i + djVi)'^dr of the resulting 
quadrupolar velocity field in the electrolyte. 

When comparing the results for the strength of the 
ICEO flow around the cylinder obtained by the analytical 
model with those obtained by a numerical solution of 
the full equation system, where the double layer is fully 
resolved, we have noted significant discrepancies. These 
discrepancies, which are described in the following, have 
become the primary motivation for the study presented 
in this paper. 

First, in the full double-layer resolving simulation we 



determined the value P* 



moch(-^o) — J|t-|=_Ro 

of the mechanical input power, where Rq is the radius 
of a cylinder surface placed co-axially with the metallic 
cylinder. Then, as expected due to the electrical forces 
acting on the net charge in the double layer, we found 
that P^cch(^o) varied substantially as long as the inte- 
gration cylinder surface was inside the double layer. For 
Rq fa Uc + 6Ad the mechanical input power stabilized at 
a certain value. However, this value is significantly lower 
than the analytical value, but the discrepancy decreased 
for decreasing values of Ad • Remarkably, even for a quite 
thin Debye layer. Ad = 0.01 Oc, the value of the full 
numerical simulation was about 40% lower than the ana- 
lytical value. Clearly, the analytical model overestimates 
the ICEO effect, and the double-layer width must be ex- 
tremely thin before the simple analytical model agrees 
well with the full model. 

A partial explanation of the quantitative failure of the 
analytical slip velocity model is the radial dependence 
of the tangential field combined with the spatial ex- 
tent of the charge density pci of the double layer. In 
the Debye-Hiickel approximation iSy and pc\ around the 
metallic cylinder of radius ac become 
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where Ki is the decaying modified Bessel function of 
order 1. The slowly varying part of is given by 
Eq [l -I- (flc/ r)^] sin 9. For very thin double layers it is well 
approximated by the r-independent expression 2Eq sin6', 
while for wider double layers, the screening charges sam- 
ple the decrease of as a function of the distance from 
the cylinder. Also tangential hydrodynamic and osmotic 
pressure gradients developing in the double layer may 



contribute to the lower ICEO strength when taking the 
finite width of the double layer into account. 

In this work we analyze quantitatively the impact 
of a finite Debye length on the kinetic energy of the 
fiow rolls generated by ICEO for three different models: 
(i) The full nonlinear electrokinetic model (FN) with a 
fully resolved double layer, (ii) the linear slip-velocity 
model (LS), where electrostatics and hydrodynamics are 
completely decoupled, and (iii) a nonlinear slip-velocity 
model (NSL) including the double layer charging through 
ohmic currents from the bulk electrolyte and the surface 
conduction in the Debye layer. The latter two models 
are only strictly valid for infinitely thin double layers, 
and we emphasize that the aim of our analysis is to de- 
termine the errors introduced by these models neglecting 
the finite width of the double layers compared to the full 
nonlinear model resolving the double layer. We do not 
seek to provide a more accurate description of the physics 
in terms of extending the modeling by adding, say, the 
Stern layer (not present in the model) or the steric effects 
of finite-sized ions (not taken into account). 



II. MODEL SYSTEM 

To keep our analysis simple, we consider a single un- 
biased metallic electrode in a uniform, external electric 
field. The electrode of width 2a and height h is placed 
at the bottom center, —a < x < a and 2: = 0, of a 
square 2L x 2L domain in the xz-plane filled with an 
electrolyte, see Fig. [T] The system is unbounded and 
translational invariant in the perpendicular y-direction. 
The uniform electric field, parallel to the surface of the 
center electrode, is provided by biasing the driving elec- 
trodes placed at the edges x — ±L with the DC volt- 
ages ±Vo, respectively. This anti-symmetry in the bias 
voltage ensures that the constant potential of the center 
electrode is zero. A double layer, or a Debye screening 
layer, is induced above the center electrode, and an ICEO 
flow is generated consisting of two counter-rotating flow 
rolls. Electric insulating walls at z = (for |a:;| > a) and 
at z = 2i confine the domain in the z-direction. The 
symmetry of the system around a; = is exploited in the 
numerical calculations. 



III. FULL NONLINEAR MODEL (FN) 

We follow the usual continuum approach to the elec- 
trokinetic modeling of the electrolytic microchamber and 
treat only steady-state problems. For simplicity we con- 
sider a symmetric, binary electrolyte, where the positive 
and negative ions with concentrations c+ and c_ , respec- 
tively, have the same diffusivity D and charge number Z. 
Using the ideal gas model for the ions, an ion is affected 
by the sum of an electrical and an osmotic force given 
by F± = TZeVcp - {k^T/c±) Vc±. Here e is the ele- 
mentary charge, T is the absolute temperature and fee is 
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densities are governed by the Nernst-Planck equation 



FIG. 1: A sketch of the square 2L x 2L electrolytic microcham- 
ber in the xz-plane. The external voltage ±Vb is applied to 
the two electrodes (thick black lines) ai x — ±L, respectively. 
It induces two counter-rotating flow rolls (curved black ar- 
rows) by electro-osmosis over the un-biased metallic center 
electrode of length 2a and height h placed at the bottom wall 
around {x,z) = (0,0). The spatial extent of the flow rolls is 
represented by the streamline plot (thin black curves) drawn 
as equidistant contours of the flow rate. The inset is a zoom-in 
on the right half, < a; < a, of the un-biased center electrode 
and the nearby streamlines. 



Boltzmann's constant. Assuming a complete force bal- 
ance between each ion and the surrounding electrolyte, 
the resulting body force density /ion = '^i=±CiFi, ap- 
pearing in the Navier-Stokes for the electrolyte due to 
the forces acting on the ions, is 



-Ze{c 



)V<j)~kBTV{c+ 



(3) 



As the second term is a gradient, namely the gradient of 
the osmotic pressure of the ions, it can in the Navier- 
Stokes equation be absorbed into the pressure gradient 
Vp = Vpdyn + Vpos, which is the gradient of the sum of 
hydrodynamic pressure and the osmotic pressure. Only 
the electric force is then kept as an explicit body force. 
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where the first term expresses ionic diffusion and the sec- 
ond term ionic electromigration due to the electrostatic 
potential (j>. The last term expresses the convective trans- 
port of ions by the fluid velocity field v. 

The electrostatic potential is determined by the charge 
density pd — Ze{c^ ^ c_) through Poisson's equation 



V ■ (£V0) = -Pel, 



(6) 



where e is the fluid permittivity, which is assumed con- 
stant. The fluid velocity field v and pressure field p are 
governed the the continuity equation and the Navier- 
Stokes equation for incompressible fluids, 

V ■v = 0, (7a) 
p„,{v ■V)v = -Vp + TjW^v - pdV<j>, (7b) 

where and 77 are the fluid mass density and viscosity, 
respectively, both assumed constant. 



B. Dimensionless form 

To simplify the numerical implementation, the gov- 
erning equations are rewritten in dimensionless form, as 
summarized in Fig. [21 using the characteristic parameters 
of the system: The geometric half-length a of the elec- 
trode, the ionic concentration cq of the bulk electrolyte, 
and the thermal voltage (po = k-QT/{Ze). The character- 
istic zeta-potential C of the center electrode, i.e. its in- 
duced voltage, is given as the voltage drop along half of 
the electrode, C, = (a/L)Vo, and we introduce the dimen- 
sionless zeta-potential a as C = a^o, or a = (aVo)/(i^o)- 
The characteristic velocity uo is chosen as the Helmholtz- 
Snioluchowski slip velocity induced by the local electric 
field E ~ (/a, and finally the pressure scale is set by the 
characteristic microfluidic pressure scale po = rjuo/a. In 
summary. 
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The new dimensionless variables (denoted by a tilde) thus 
become 



A. Bulk equations 

Neglecting bulk reactions in the electrolyte, the ionic 
transport is governed by the particle conservation 



V ■ J± = 0, 



(4) 



where J± is the flux density of the two ionic species. As- 
suming the electrolytic solution to be dilute, the ion flux 
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Ci = — , 
Co 



P=^. (9) 
uo Po 



To exploit the symmetry of the system and the binary 
electrolyte, the governing equations are reformulated in 
terms of the average ion concentration c = (c+ -|- c_)/2 
and half the charge density p = (c+ — c_)/2. Correspond- 
ingly, the average ion flux density Jc = {J+ + J-)/'2, and 
half the current density Jp = {J^ — JJ)/2 are introduced. 
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dzC = 0, Vx = Vz = 
dzp = 0, 9,0 = 



dzC = 0, 
dzP = 0, 



Vx = Vz 

= 



-L 








ReVjdjVt = djGij - -^pdicf) 



Pe VjdjC 



dj{djC + pdjcj) 
dj {djp + cdj (p) = Pe Vj dj p 



dzC = 0, 
dzP = 0, 



Vx = Vz = 
dz<f>=0 







1, 

0, 







Vz=0 



FIG. 2: The governing equations (without box) and boundary conditions (with boxes, arrows points to specific boundaries) in 
dimensionless form (the tiide is omitted for clarity) for the entire quadratic 2L x 2L domain (not shown in correct aspect ratio) 
bisected into two symmetric halves. Only the right half {x > 0) of the domain is included in the simulations. The boundaries 
are the surface of the un- biased center electrode (black rectangle), the solid insulating walls (dark gray lines), the external 
electrode (black line), and the symmetry line (dashed black line). 



Thus, the resuhing full system of coupled nonlinear equa- 
tions takes the following form for the ionic fields 



V • = V • Jp = 0, 



(10a) 

Jc = -pV(j) - Vc + Pecv, (10b) 
Jp = ~cV4>-'Vp+ Pepv, (10c) 

Pe-^, (lOd) 



while the electric potential obeys 
V^=-^p, 

and finally the fiuid fields satisfy 
V • f; = 0, 



(11) 



(12a) 



Re(v ■V)v = -Vp + V^v - ^^V(t), (12b) 



Re = 



puoa 



(12c) 



Here the small dimensionless parameter e = Ad/o has 
been introduced, where Ad is the Debye length. 



Ad 1 / skj^T 



a aV 2(Ze)-^co 



(13) 



We note that the dimensionless form of the osmotic force, 
the second term in Eq. (O, is f°^^-^ = -(l/e2a^)Vc. 



C. Boundary conditions 

We exploit the symmetry around a; = and consider 
only the right half (0 < a; < L) of the domain, see Fig. [21 
As boundary conditions on the driving electrode we take 



both ion concentrations to be constant and equal to the 
bulk charge neutral concentration. Correspondingly, the 
charge density is set to zero. Consequently, we ignore 
all dynamics taking place on the driving electrode and 
simply treat it as an equipotential surface with the value 
Vq- We set a no-slip condition for the fluid velocity, and 
thus at X = i we have 

£=1, p = 0, 0=-^=a-, v = 0. (14) 
00 a 

On the symmetry axis {x = 0) the potential and the 
charge density must be zero due to the anti-symmetry of 
the applied potential. Moreover, there is neither a fluid 
flux nor a net ion flux in the normal direction and the 
shear stresses vanish. So at x = we have 

= 0, n • Jc 0, p = 0, (15a) 
t ■ a- ■ n = 0, fi ■ V = 0, (15b) 

where the stress tensor is (cr)ifc = —pSik + r]{diUk + dkUi) , 
and n and t are the normal and tangential unit vectors, 
respectively, which in 2D, contrary to 3D, are uniquely 
defined. The constant potential on the un-biased metallic 
electrode is zero due to symmetry, and on the electrode 
surface we apply a no-slip condition on the fluid velocity 
and no-current condition in the normal direction. So on 
the electrode surface we have 



n • Jc = 0, n • Jp = 0, 



= 0, v = 0. 



(16) 



On the solid, insulating walls there are no fluxes in the 
normal direction, the normal component of the electric 
fleld vanishes and there are no-slip on the fluid velocity. 



n • Jc = 0, ri • Jo = 0, n ■ V0 = 0, v = 0. 



(17) 



A complete overview of the governing equations and 
boundary conditions is given in Fig. [21 
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dzC=0, Vx = Vz = 0] 

d,p = 0, 5,0 = 



Vx = 0, dxVz = 
= 0, dxc = 
p^O, dxP + 2dx4> = 



d^c = 0, Vx = v^^Q 
d^p + 29,0 = 0, = 



-L 







djVj = 

ReVjdjUi = djaij - ■^e"^'^ smh{^) dicj) 
^20= -^e«/2sinh(|) 
-9|c + Pev,d,c = i [(9,c)2 + id.pf] + (9,0)(9,p) 



+ 20) = (9,c)(9,p) + (9,0)(9,c) - PeT;,9,p 



a^c = 0, -y^ = = 
= 0, a,0 = 



C = 0, Ui: 

p = 0, 



FIG. 3: The governing equations (without box) and boundary conditions (with boxes) in dimensionless form (the tilde is 
omitted) using the logarithmic concentrations (denoted by a breve) of Eq. (|18p . Otherwise the figure is identical to Fig. [J] 



D. The strongly nonlinear regime 



IV. SLIP- VELOCITY MODELS 



At high values of the induced ^-potential, the concen- 
trations of counter- and co-ions acquire very large and 
very small values, respectively, near the center electrode. 
Numerically this is problematic. The concentration ratio 
becomes extremely large and the vanishingly small con- 
centration of co-ions is comparable to the round-off error 
and may even become negative. However, this numerical 
problem can be circumvented by working with the loga- 
rithms (marked by a breve accent) of the concentration 
fields, c± = log(c±/co). By inserting 



c± = Co exp ( c± I 



(18) 



in the governing equations (15|) , ([6]) , and (j7bp , a new equiv- 
alent set of governing equations is derived. The symme- 
try is exploited by defining the symmetric c = c+ + C- 
and antisymmetric p = c+ — c_ combination of the log- 
arithmic fields and the corresponding formulation of the 
governing equations is 



-20; 



Pe V ■ Vc — 



(Vc)2 + (Vp)2 



Pev -Vp-Vc-Vp-Vc, 



1 



- V0 • Vp, 
(19a) 
•Vp, (19b) 

(19c) 



Re(v ■ V)v = -Vp 



1 



^,^,e^/2sinh(^|V, 



(19d) 



while the continuity equation remains the same as in 
Eq. (|12ap . The governing equations and boundary con- 
ditions for the logarithmic fields (breve-notation) is sum- 
marized in Fig. [3] This transformation serves to help 
linearize solutions of the dependent variables, c and p, 
at the expense of introducing more nonlinearity into the 
governing equations. 



The numerical calculation of ICEO flows in microflu- 
idic systems is generally connected with computational 
limitations due to the large difference of the inherent 
length scales. Typically, the Debye length is much 
smaller than the geometric length scale. Ad ^ a, mak- 
ing it difflcult to resolve both the dynamics of the Debye 
layer and the entire microscale geometry with the avail- 
able computer capacity. Therefore, it is customary to 
use slip-velocity models, where it is assumed that the 
electrodes are screened completely by the Debye layer 
leaving the bulk electrolyte charge neutral. The dynam- 
ics of the Debye layer is modeled separately and applied 
to the bulk fluid velocity through an effective, so-called 
Helmholtz-Smoluchowski slip velocity condition at the 
electrode surface. 



■Whs 



V 



(20) 



where C is the zeta potential at the electrode surface, 
and is the electric fleld parallel to the surface. Re- 
gardless of the modeled dynamics in the double layer the 
slip- velocity models are only strictly valid in the limit of 
infinitely thin double layers Ad -C a. 



A. The linear slip-velocity model (LS) 

The double-layer screening of the electrodes leaves the 
bulk electrolyte charge neutral, and hence the governing 
equations only include the potential 0, the pressure field 
p and the flow velocity fleld v. In dimensionless form 
they become. 



V^0 = O, 
Re(v ■ 'V)v = -Vp + V^v, 
V-v = 0. 



(21a) 
(21b) 
(21c) 



6 



The electrostatic problem is solved independently of the 
hydrodynamics, and the potential is used to calculate 
the effective slip velocity applied to the fluid at the un- 
biased electrode surface. The boundary conditions of the 
potential and fluid velocity are equivalent to the condi- 
tions applied to the full non-linear system, except at the 
surface of the un-biased electrode. Here, the normal com- 
ponent of the electric fleld vanishes, and the effective slip 
velocity of the fluid is calculated from the electrostatic 
potential using ( = —(/> and E^^ = — [(t • V)0] t, 



where a is the bulk 3D conductivity and 



1 



■Whs 



(t-V) 



(22a) 
(22b) 



This represents the simplest possible, so-called linear 
slip-velocity model; a model which is widely applied as 
a starting point for numerical simulations of actual mi- 
crofluidic systems [l^, [2l[ . In this simple model all the 
dynamics of the double layer has been neglected, an 
assumption known to be problematic when the voltage 
across the electrode exceeds the thermal voltage. 



B. The nonlinear slip-velocity model (NLS) 

The linear slip- velocity model can be improved by tak- 
ing into account the nonlinear charge dynamics of the 
double layer. This is done in the so-called nonlinear slip- 
velocity model, where, although still treated as being in- 
finitely thin, the double layer has a non-trivial charge 
dynamics with currents from the bulk in the normal di- 
rection and currents flowing tangential to the electrode 
inside the double layer. For simplicity we assume in the 
present nonlinear model that the neutral salt concentra- 
tion Co is uniform. This assumption breaks down at high 
zeta potentials, where surface transport of ionic species 
can set up gradients in the salt concentrations leading 
to chemi-osmotic flow. In future more complete studies 
of double layer charge dynamics these effects should be 
included. 

The charging of the double layer by the ohmic 
bulk current is assumed to happen in quasi-equilibrium 
characterized by a nonlinear differential capacitance 
Cdi given by the Gouy-Chapmann model, Cdi — 
e cosh[zeC/(2fcBr)]/AD, which in the the low-voltage, lin- 
ear Debye-Hiickel regime reduces to Cdi = s/Xu- Ignor- 
ing the Stern layer, the zeta-potential is directly propor- 
tional to the bulk potential right outside the double layer, 

C = -0. 

The current along the electrode inside the Debye layer 
is described by a 2D surface conductance CTs, which for a 
binary, symmetric electrolyte is given by pj] 



as = 4 Ad 17(1 + m) sinh" 



(ZeC_\ 



(23) 



m 



kBT 
't]D \ Ze 



(24) 



is a dimensionless parameter indicating the relative con- 
tribution of electroosmosis to surface conduction. In 
steady state the conservation of charge then yields [11] 



= n • {aVc, 



(25) 



where the operator = t{t ■ V) is the gradient in the 
tangential direction of the surface. 

Given the length scale a of the electrode, the strength 
of the surface conductance can by characterized by the 
dimensionless Dukhin number Du deflned by 



Du 



(Js _ 4Ad 
acr a 



to) sinh 



(26) 



Conservation of charge then takes the dimensionless form 
Q = n-{Vij}) + Vs-[DuVs-$\, (27) 

and this effective boundary condition for the potential 
on the flat electrode constitutes a ID partial differen- 
tial equation and as such needs accompanying boundary 
conditions. As a boundary condition the surface flux is 
assumed to be zero at the edges of the electrode, 



{t-V) 



I a;— ±a 



0, 



(28) 



which is well suited for the weak formulation we employ 
in our numerical simulation as seen in Eq. (|34p . 



V. NUMERICS IN COMSOL 

The numerical calculations are performed using the 
commercial finite-element-method software COMSOL 
with second-order Lagrange elements for all the fields ex- 
cept the pressure, for which first-order elements suffices. 
We have applied the so-called weak formulation mainly 
to be able to control the coupling between the bulk equa- 
tions and the boundary constraints, such as Eqs. (|22b|) 
and (|25p , in the implementation of the slip- velocity mod- 
els in script form. 

The Helmholtz-Smoluchowski slip condition poses a 
numerical challenge because it is a Dirichlet condition 
including not one, but up to three variables, for which 
we want a one-directional coupling from the electrostatic 
field (t> to the hydrodynamic fields v and p. We use the 
weak formulation to unambiguously enforce the bound- 
ary condition with the explicit introduction of the re- 
quired hydrodynamic reaction force / on the un-biased 
electrode 



f =CT h. 



(29) 



The X and z components of Navier-Stokes equation are 
multiplied with test functions and Uz, respectively. 
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and subsequently integrated over the whole domain Q. 
Partial integration is then used to move the stress tensor 
contribution to the boundaries dfl, 



0=1 UiCTijUjAs 

dn 



[{djUt)aij + UiBi] da, (30) 



where Bi = Re (vjdj) Vi + p{di<j}) / {e a ) . The boundary 
integral on the un-biased electrode dVluc is rewritten as 



ds = 



[uifi + Fi{vi - V}is.i)]ds, (31) 



where Fi are the test functions belonging to the compo- 
nents fi of the reaction force /. These test functions are 
used to enforce the Helmholtz-Smoluchowski slip bound- 
ary condition consistently. This formulation is used for 
both slip- velocity models. 

In the nonlinear slip-velocity model the Laplace equa- 
tion (|21ap is multiplied with the electrostatic test func- 
tion $ and partially integrated to get a boundary term 
and a bulk term 







$ (9,0) n,ds - / (5,$) id,(f>) da. (32) 
on Jn 



The boundary integration term on the electrode is sim- 
plified by substitution of Eq. ([25ll which results in 



$ {d,(t)) n,ds = $ [Ud, {Duijdj<j))] ds. 

.(33) 

Again, the resulting boundary integral is partially inte- 
grated, which gives us explicit access to the end-points 
of the un-biased electrode. This is necessary for applying 
the boundary conditions on this ID electrode, 



/ $ Ud,{Duijdj(l)) ds 



x—-\-a 



{Ud^^)Du {ijdj(p) ds 



(34) 



The no-flux boundary condition can be explicitly in- 
cluded with this formulation. Note that in both slip- 
velocity models the zeta-potential is given by the poten- 
tial just outside the Debye layer, ( = —0, and it is there- 
fore not necessary to include it as a separate variable. 

The accuracy and the mesh dependence of the simu- 
lation as been investigated as follows. The comparison 
between the three models quantifies relative differences 
of orders down to 10^^, and the convergence of the nu- 
merical results is ensured in the following way. COMSOL 
has a build-in adaptive mesh generation technique that 
is able to refine a given mesh so as to minimize the error 
in the solution. The adaptive mesh generator increases 
the mesh density in the immediate region around the 
electrode to capture the dynamics of the ICEO in the 
most optimal way under the constraint of a maximum 
number of degrees of freedom (DOFs). For a given set 



of physical parameters, the problem is solved each time 
increasing the number of DOFs and comparing consec- 
utive solutions. As a convergence criterium we demand 
that the standard deviation of the kinetic energy relative 
to the mean value should be less than a given threshold 
value typically chosen to be around 10~^. All of the sim- 
ulations ended with more than 10^ DOFs, and the ICEO 
flow is therefore sufficiently resolved even for the thinnest 
double layers in our study for which e = 10^^. 



VI. RESULTS 

Our comparison of the three numerical models is pri- 
marily focused on variations of the three dimensionless 
parameters e, a, and (3 relating to the Debye length Ad, 
the applied voltage Vb, and the height h of the electrode, 
respectively, 



Ad 
a 



aVo 
L4>o 



h 

a 



(35) 



As mentioned in Sec. HI the strength of the generated 
ICEO flow can be measured as the mechanical power in- 
put Pmcch exerted on the electrolyte by the slip-velocity 
just outside the Debye layer or equivalently by the ki- 
netic energy dissipation Pkin m the bulk of the electrolyte. 
However, both these methods suffers from numerical in- 
accuracies due to the dependence of both the position 
of the integration path and of the less accurately deter- 
mined velocity gradients in the stress tensor er. To ob- 
tain a numerically more stable and accurate measure, we 
have chosen in the following analysis to characterize the 
strength of the ICEO flow by the kinetic energy Skin of 
the induced flow fleld 



kin 



(36) 



which depends on the velocity fleld and not its gradi- 
ents, and which furthermore is a bulk integral of good 
numerical stability. 



A. Zero height of the un-biased center electrode 

We assume the height h of the un-biased center elec- 
trode to be zero, i.e. /? = 0, while varying the Debye 
length and the applied voltage through the parameters 
e and a. We note that the linear slip-velocity model 
Eqs. ([2T]l and ((22|) is independent of the dimensionless 
Debye length e. It is therefore natural to use the kinetic 
energy E\^f^ of this model as a normalization factor. 

In the lin-log plot of Fig. [1] we show the kinetic energy 
E^^^ and E™ normalized by E^f^ as a function of the 
inverse Debye length 1/e for three different values of the 
applied voltage, a = 0.05, 0.5 and 5, ranging from the 
linear to the strongly nonlinear voltage regime. 
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10° 10^ 10^ 10^ 10^ 

1/e 

FIG. 4: The total induced kinetic energy -Ekin^ (gray dashed) 
and i?kin (black) for the nonlinear slip- velocity model and the 
full model, respectively, relative to (horizontal black line) 
of the linear slip- velocity model as a function of dimensionless 
inverse Debye length 1/e. Each are shown for three values of 
the dimensionless applied voltage a — 0.05, 0.5 and 5. The 
value of e decreases from 1 to lO"** going from left to right. 



We first note that in the limit of vanishing Debye 
length (to the right in the graph) all models converge 
towards the same value for all values of the applied volt- 
age a. For small values of a the advanced slip-velocity 
model E^^^ is fairly close to the linear slip- velocity model 
-^kin' t>ut as a increases, it requires smaller and smaller 
values of e to obtain the same results in the two models. 
In the linear regime a = 0.05 a deviation less than 5% 
is obtained already for e < 1. In the nonlinear regime 
a = 0.5 the same deviation requires e < 10~^, while 
in the strongly nonlinear regime e < 10""* is needed to 
obtain a deviation lower than 5%. 

In contrast, it is noted how the more realistic full model 
E™ deviates strongly from E^^ for most of the displayed 
values of e and a. To obtain a relative deviation less 
than 5% in the linear (a — 0.05) and nonlinear {a — 0.5) 
regimes, a minute Debye length of e < lO^'^ is required, 
and in the strongly nonlinear regime the 5% level it not 
reached at all. 

The deviations are surprisingly large. The Debye 
length in typical electrokinetic experiments is Ad — 
30 nm. For a value of e = 0.01 this corresponds to 
an electrode of width 2x3 unr — 6 pm, comparable 

to those used in Refs. Bin p. In Fig. a we see that 

for a = 5, corresponding to a moderate voltage drop of 
0.26 V across the electrode, the linear slip- velocity model 
overestimates the ICEO strength by a factor 1/0.4 — 2.5. 
The nonlinear slip-model does a better job. For the same 
parameters it only overestimates the ICEO strength by 
a factor 0.5/0.4 = 1.2. 

For more detailed comparisons between the three mod- 




10° 10^ 10^ 10^ lO'' 

FIG. 5: The difference between the induced kinetic energies 
-Ejin and -Ekin^ of the linear and nonlinear slip- velocity mod- 
els, respectively, relative to the full model -E^in a function 
of the inverse Debye length 1/e. for three different applied 
voltages a — 0.05, 0.5, 5. 




10-1 10° 



a 

FIG. 6: The difference between the induced kinetic energies 
and -E^in^ of the linear and nonlinear slip-velocity mod- 
els, respectively, relative to the full model E™ as a function of 
the voltage bias a for three different Debye layer thicknesses 
e = 1.8 X 10"^ 10"^ 10"\ 



els the data of Fig. H] are plotted in a different way 
in Fig. [5] Here the overestimates (E^^^/ E™) - 1 and 
(i?^[;'^/i?™) — 1 of the two slip- velocity models relative 
to the more correct full model are plotted in a log-log 
plot as a function of the inverse Debye length 1/e for 
three different values of the applied voltage. It is clearly 
seen how the relative deviation decreases proportional to 
e as e approaches zero. 

Finally, in Fig.glthe relative deviations {E^^^jEl^)-l 
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and (-Ej^Ji,'^ / E™) — 1 are plotted versus the voltage a in a 
log-log plot. For any value of the applied voltage a, both 
slip- velocity models overestimates by more than 100% for 
large Debye lengths e = lO"'^ and by more than 10% for 
e = 10"^. For the minute Debye length Ad = 1-8 x lO^'^ 
the overestimates are about 3% in the linear and weakly 
nonlinear regime a < 1, however, as we enter the strongly 
nonlinear regime with a = 5 the overestimation increases 
to a level above 10%. 



and /? = 0.001 are close. As the height is increased to 
P — 10^^ we note that the strength of the ICEO is in- 
creased by 20%— 25% as f3 > e. This tendency is even 
stronger pronounced for the higher electrode f3 = 10~^. 
Here the ICEO strength is increased by approximately 
400% for a large range of Debye lengths. 



C. Thermodynamic efficiency of the ICEO system 



B. Finite height of the un-biased electrode 

Compared to the full numerical model, the slip- velocity 
models are convenient to use, but even for small Debye 
lengths, say Ad = O.Ola, they are prone to significant 
quantitative errors as shown above. Similarly, it is of 
relevance to study how the height of the un-biased elec- 
trode influences the strength of the ICEO flow rolls. In 
experiments the thinnest electrodes are made by evapo- 
ration techniques. The resulting electrode heights are of 
the order 50 nm — 200 nm, which relative to the typi- 
cal electrode widths a « 5 pm results in dimensionless 
heights 10-3 < P < 10"^ 

In Fig. [7] is shown the results for the numerical calcu- 
lation of the kinetic energy £^™(e,/3) using the full nu- 
merical model. The dependence on the kinetic energy of 
the dimensionless Debye length e — Xd/q and the dimen- 
sionless electrode height (3 — h/ a is measured relative to 
the value £^™(e,/9) of the infinitely small Debye length 
for an electrode of zero height. For small values of the 
height no major deviations are seen. The curve for /3 = 



Z a 
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FIG. 7: The difference between the induced kinetic en- 
ergies E^^{t,(3) of the full model at finite Debye length 
and electrode height relative to the full model i?^n(0,0) at 
zero Debye length and zero electrode height as a function 
of the inverse Debye length 1/e for four electrode heights 
/3 = 0,10-^10-^10-^ 



Conventional electro-osmosis is known to have a low 
thermodynamic efficiency defined as the delivered me- 
chanical pumping power relative to the total power de- 
livered by the driving voltage. Typical efhciencies are of 
the order of 1% . 26], while in special cases an efhciency of 
5.6% have been reported |23] . In the following we provide 
estimates and numerical calculations of the correspond- 
ing thermodynamic efficiency of the ICEO system. 

The applied voltage drop 2Vb = 2EqL across the sys- 
tem in the x-direction is written as the average electrical 
field Eo times the length 2L, while the electrical cur- 
rent is given by / = WHaEo, where W and H is the 
width and height in the y- and z-direction, respectively, 
and a = De/X^ = e/tu is the conductivity written in 
terms of the Debye time td = Aq/Z?. The total power 
consumption to run the ICEO system is thus 



Pt, 



2Vo X / 



4 

TD 



LWH. 



(37) 



This expression can be interpreted as the total energy, 
■iei?o X LWH, stored in the average electrical field of 
the system with volume LWH multiplied by the charac- 
teristic electrokinetic rate 4/td. 

The velocity-gradient part of the hydrodynamic stress 



tensor is denoted cr, i.e. (cr)^ 



-qidiVj + djVi). In 



terms of cr, the kinetic energy dissipation Pkin neces- 
sary to sustain the steady-state flow rolls is given by 
-Pkin = ^ /o^ Tr(<T^). In the following esti- 

mate we work in the Debye-Hiickel limit for an electrode 
of length 2a, where the induced zeta potential is given 
by Cind — aEo and the radius of each flow roll is approx- 
imately a. In this limit the electro-osmotic slip velocity 
Weo and the typical size of the velocity gradient \diVj \ are 



En = 



ea 



E, 



\diVj\ w — ^ - Eq. 
a J] 



(38a) 
(38b) 



Thus, since the typical area covered by each flow roll is 
Tra^ , we obtain the following estimate of Pkin, 



W , 



V- 
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sE^ (1 



{^eEo) ira^W. (39) 



Here the power dissipation can be interpreted as the en- 
ergy of the electrical field in the volume na^W occupied 
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by each flow roll multiplied by an ICEO rate given by the 
electric energy density eEq divided by the rate of viscous 
energy dissipation per volume given by rj. 

The thermodynamic efficiency can now be calculated 
as the ratio Pkin/^tot using Eqs. ((37)) and p9|) . 



Pi 



kin 



2.4 X 10" 



(40) 



This efficiency is the product of the ratio between the 
volumes of the flow rolls and the entire volume multiplied 
and the ratio of the electric energy density in the viscous 
energy density 77/10. The value is found using L ^ H ^ 
15a — 0.15 mm, Eq = 2.5 kV/m, and Ad = 20 nm, 
which is in agreement with the conventional efficiencies 
for conventional electro-osmotic systems quoted above. 



VII. CONCLUSION 

We have shown that the ICEO velocities calculated 
using the simple zero-width models significantly overes- 
timates those calculated in more realistic models taking 
the finite size of the Debye screening length into ac- 
count. This may provide a partial explanation of the 
observed quantitative discrepancy between observed and 
calculated ICEO velocities. The discrepancy increases 
substantially for increasing e, i.e. in nanofluidic systems. 

Even larger deviations of the ICEO strength is calcu- 
lated in the full numerical model when a small, but finite 
height of the un-biased electrode is taken into account. 

A partial explanation of the quantitative failure of the 
analytical slip velocity model is the decrease of the tan- 
gential electric field as a function of the distance to the 



surface of the polarized ICEO object combined with the 
spatial extent of the charge density of the double layer. 
Also tangential hydrodynamic and osmotic pressure gra- 
dients developing inside the double layer may contribute 
to the lowering ICEO strength when taking the finite 
width of the double layer into account. The latter may 
be related to the modification of the classical Helmholtz- 
Smoluchowski expression of the slip- velocity obtained by 
adding a term proportional to the gradient of the salt 
concentration c [28] . 

Our work shows that for systems with a small, but 
non-zero Debye length of 0.001 to 0.01 times the size of 
the electrode, and even when the Debye-Hiickel approx- 
imation is valid, a poor quantitative agreement between 
experiments and model calculations must be expected 
when applying the linear slip-velocity model based on 
a zero Debye-length. It is advised to employ the full 
numerical model of ICEO, when comparing simulations 
with experiments. 
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